\documentclass[aps,twocolumn]{revtex4}
\usepackage{graphicx}
\usepackage{amssymb,amsfonts,amsmath}
\usepackage{chemarr}
\usepackage{bm}
\usepackage{pslatex}
%\usepackage{mathptmx}
\usepackage{xfrac}

\newcommand{\mymat}[1]{\boldsymbol{#1}}
\newcommand{\mytrn}[1]{{\!\!~^{\mathsf{t}}{#1}}}
\newcommand{\myvec}[1]{\overrightarrow{#1}}
\newcommand{\half}{\sfrac{1}{2}}
\newcommand{\onethird}{\sfrac{1}{3}}
\newcommand{\twothirds}{\sfrac{2}{3}}
\newcommand{\q}{\vec{q}}
\newcommand{\dq}{\dot{\q}}
\newcommand{\ddq}{\ddot{\q}}
\newcommand{\C}{\vec{C}}
\newcommand{\J}{\mymat{J}}
\newcommand{\dJ}{\dot{\J}}
\newcommand{\tJ}{\mytrn{\J}}
\newcommand{\G}{\vec{G}}
\newcommand{\W}{\mymat{W}}
\newcommand{\A}{\mymat{A}}
\newcommand{\JW}{\hat{\mymat{J}}}
\newcommand{\tJW}{\mytrn{\JW}}

\begin{document}

\title{Some particles dynamics, with constraints}
\maketitle
\section{Setup}
Let $N$ be the number of particles with coordinates $\q$.	
we assume that those particles have a set of external forces $\vec{F}$.
We also assume that there is a set of $M$ holonomic, stationary constraints
\begin{equation}
	\C\left(\q\right) = \vec{0}
\end{equation}
which define the legal coordinates.
We also define the legal velocities by
\begin{equation}
	\J \dq = \vec{0}.
\end{equation}
We assume that those internal constraints yields internal forces which are
a combination of the Jacobian lines (or of the rows of the Jacobian transpose).
Let us devise the velocity Verlet algorithm for this system, during a time step $\tau$.

\section{Algorithm}

\subsection{Initial phase space and forces}
We avec $\q_n$ and $\dq_n$ such that
\begin{equation}
	\C(\q_n) = \vec{0}
\end{equation}
and
\begin{equation}
	\J_n \dq_n = \vec{0}.
\end{equation}
We also compute
\begin{equation}
	\vec{F}_n = \vec{F}\left(\q_n,\dq_n\right) 
\end{equation}

\subsection{Predicted Position}
The predicted position is
\begin{equation}
	\q_{n+\half} = \q_{n} + \tau \dq_n + \dfrac{1}{2} \tau^2 \W \vec{F}_n
\end{equation}

\subsection{Corrected Position}
We have to come back to the legal position using a set of Lagrange multipliers so that
\begin{equation}
	\C\left( \q_{n+\half} + \W\tJ_n \lambda_n \right) = \vec{0}
\end{equation}
We use a PseudoNewton approach using
$$
	\C\left( \q_k + \W\tJ_n \lambda_k \right) \simeq \C\left( \q_k \right) + \J_n\W\tJ_n \lambda_k.
$$
We define $\A_n = \J_n\W\tJ_n$ which must be positive and symmetric, otherwise the phase space is singular.\\

We end up with
\begin{equation}
	\q_{n+1} = \q_{n+\half} + \W\tJ_n \lambda_n = \q_{n+\half} + \delta\q_n.
\end{equation}
An estimate of the internal forces is
$$
	\vec{G}_n = \dfrac{2}{\tau^2} \delta\q_n
$$

\subsection{Predicted half-velocities}
The first contribution to the velocities is
\begin{equation}
	\dq_{n+\sfrac{1}{4}} = \dq_n + \frac{1}{2}\tau\W \vec{F}_n+ \dfrac{1}{\tau}\delta\q_n
\end{equation}

\subsection{Half Velocities Legalisations}
We can compute $J_{n+1}$.
\begin{equation}
	\dq_{n+\half} = \dq_{n+\sfrac{1}{4}} + \W\tJ_{n+1} \vec{\mu}_n
\end{equation}
and
\begin{equation}
	\vec{0} = \J_{n+1}\dq_{n+\sfrac{1}{4}} + \underbrace{\J_{n+1}\W\tJ_{n+1}}_{A_{n+1}} \vec{\mu}_n
\end{equation}

\subsection{Velocities Correction}
We can now compute
\begin{equation}
	\vec{F}_{n+\half} = \vec{F}(\q_{n+1},\dq_{n+\half})
\end{equation}
and compute
\begin{equation}
	\dq_{n+\sfrac{3}{4}} = \dq_{n+\half} + \dfrac{1}{2}\tau \W\vec{F}_{n+\half}
\end{equation}
\subsection{Final Velocities Correction}
We do the same than for the half-velocities correction.

\end{document}


\end{document}

